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We investigate the instability of the saturated ferromag- 
netic ground state (Nagaoka state) in the Hubbard model on 
various lattices in dimensions d = 2 and d = 3. A variational 
resolvent approach is developed for the Nagaoka instability 
both for U = oo and for U < oo which can easily be evalu- 
ated in the thermodynamic limit on all common lattices. Our 
results significantly improve former variational bounds for a 
possible Nagaoka regime in the ground state phase diagram of 
the Hubbard model. We show that a pronounced particle-hole 
asymmetry in the density of states and a diverging density of 
states at the lower band edge are the most important features 
in order to stabilize Nagaoka ferromagnetism, particularly in 
the low density limit. 
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I. INTRODUCTION 

It is by now an often repeated fact that the so-called 
(single-band) Hubbard model was originally introduced 
to explain ferromagnetismEl u. In what followed, how- 
ever, it turned out to be rather a generic model for anti- 
ferromagnetism. Ferromagnetism seemed to require ad- 
ditional ingredients, for instance the existence of degen- 
erate bands which favor ferromagnetism based on Hund's 
rule or in the insulating case certain additional ferromag- 
netic couplings and/or correlated hopping terms. Both, 
scenarios were proven rigorously in recent years (seeu 
and references therein for the former and3 and references 
therein for the latter). 

The Hubbard model and its possible ferromagnetic 
ground state are of renewed interestafl. There are many 
works in the field based on quasi one dimensional (d = 1) 
systems triggered by the prediction of ferromagnetism 
in double minima systems at low particle densitya and 
by the numerous possibilities of analytical and numerical 
calculationsBii3 in d = 1. Exact [Galailations are possi- 
ble in infinite dimensions (d = oo)I13'Ej. For intermediate 
dimensions (1 < d <_pol numerical and approximate 
methods are employedt3il3. 

An important mile stone in the research of ferropiatt- 
netism in Hubbard models is the work of NagaokaMllia. 
It showed that at infinite local repulsion a single elec- 
tron above half-filling favors the saturated ferromagnetic 
ground state (henceforth: Nagaoka state) if the under- 
lying lattice has loops which allow interference. For bi- 
partite lattices particle-hole symmetry extends these re- 
sults to hole doping. This result reveals the beauty and 
the difficulty of the question for which lattices and for 



which fillings the Nagaoka state is the ground state. At 
U = oo, T = there is only the hopping left as a global 
energy scale. Thus there is no expansion parameter, no 
adiabatic limit, and no competition of energy scales. The 
issue is solely a question of the lattice structure, i.e. the 
possible paths on the lattice, and of the filling. 

Unfortunately, there are no extensions of Nagaoka's re- 
sult to macroscopic dopings. Qjjfc non-macroscopic num- 
bers of holes could be treatedErEJ. Therefore, we choose 
another route in the present work and investigate the sta- 
bility of the Nagaoka state towards a single spin flip. If 
such a flip lowers the energy then the Nagaoka state is 
not the ground state. Otherwise it is locally stable. The 
drawback that we treat only local stability in this way is 
not very serious. There is no indication that the tran- 
sition away from saturation should not be continuous at 
T = 0, see e.g.B 

A more serious drawback is the fact that even the single 
spin flip is too difficult a problem to be solved completely 
on finite dimensional lattices. In the limit of. infinite di- 
mensional lattices, however, it was solvecti Thereby 
it was shown that relatively simple variational ansatzes 
provide already a qualitative insight in the tendency of a 
certain lattice to have a Nagaoka state as ground state. 
Wurth et al. showed that only an extremely sophisticated 
variational ansatza yields a further reduction of the re- 
gion of possible Nagaoka state stability in comparison to 
simpler ansatzesEj. 

It is the aim of the present paper to [-extend previ- 
ous work on variational ansatzes decisively^, both in the 
completeness of the ansatzes and in the types of the lat- 
tices considered. So far, variational ansatzes considered a 
finite vicinity of the flipped spin and treated a finite num- 
ber of parameters leading to matrix eigenvalue problems. 
Here we will show that a resolvent approach is capable to 
deal implicitly with an infinite number of variational pa- 
rameters. No explicit knowledge of the variational wave 
function is required. A similar approach was used re- 
cently by OkabeO for the square lattice and the simple 
cubic lattice, too. In his work, however, the reduction of 
the resolvent to simple integrals over the density of states 
(DOS), which we succeeded to achieve in most cases, is 
lacking. 

We will present elegant simple expressions for the Na- 
gaoka instability line U CT (n) which apply to most common 
lattices. These results make it possible for everyone to 
check easily whether or not one can expect a ferromag- 
netic ground state for a given lattice. We will show that 
two main features favor the occurrence of a saturated 
ferromagnetic ground state 
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1. a highly asymmetric density of states with large 
values at the lower band edge (after particle-hole 
transformation) . 

2. non-bipartiteness of the lattice, i.e. frustration due 
to loops of three sites. 

Of course, the two points are intimately related. 

The setup of our article is as follows. In the rest of the 
Introduction we will present certain variational ansatzes 
used so far to investigate the Nagaoka state stability. In 
the subsequent section || we develop the resolvent ap- 
proach which yields simple formulae for the stability lines 
on homogeneous, isotropic lattices with nearest neighbor 
hopping. In sect. [II we present our results for various lat- 
tices in dimensions d = 2 and d — 3, namely the square, 
the simple cubic, the bee, the honeycomb, the triangular, 
the kagome, and the fee (hep) lattice. For the t-t' Hub- 
bard model on the square lattice a perturbative approach 
for small \t\ is employed as well. Sect. IV contains a sum- 
mary and a final discussion of the lattice dependence of 
saturated ferromagnetism in the Hubbard model. The 
appendices contain technical details in the derivation for 
the various lattices. 



A. Preliminary approaches 

We consider the conventional single band Hubbard 
model 

H = -t^2 a t a l° + U Y1 a ti a *J a ti a il (!) 
(bi) a i 

and calculate the spin flip energy 

Ae = {^\H-E u \^>)/{^) (2) 

where Ej\f is the energy of the Nagaoka state and is 
a variational wave function. Whenever Ae < the Na- 
gaoka state is definitely not the ground state due to the 
variational nature of our approach. At U = oo, the zero 
of Ae 00 (S) :— Ae(U — oo,5) gives the critical hole den- 
sity <5 cr above which the Nagaoka state is unstable. For 
finite U, Ae(U, 5) = leads to the Nagaoka instability 
line U CT (S) which separates a region of guaranteed insta- 
bility of the Nagaoka state (U < U CI (S)) in the phase 
diagram of the Hubbard model from a region of possible 
stability of the Nagaoka state (U > U cr (S)). In the phase 
diagrams displayed in this paper we will always represent 
the on-site repulsion U in terms of t/ rc d = U/(U + Ubr) 
where Uus = — 16e° denotes the Brinkman-Rice critical 
couplingEII. e° is the energy per particle of the saturated 
ferromagnetic state for the quarter-filled band and de- 
pends on the underlying lattice. This representation is 
chosen to render comparisons between different lattices 

A common starting point&B'BEl is defined by the 
ansatz 



|*) : = |A|~ 1/2 ]T exp(ifc b i) [a iT (a+ + / ■ £ , 



(id) 



+ .9-a+a iT ]a+|AO 



(3) 



For / = this is the Gutzwiller single spin flip wave 
function (Gw). The parameter g controls the probability 
of double occupancy. The system size is denoted by |A|. 
We use the operators a (a + ) for site diagonal fermion an- 
nihilation (creation) and c(c + ) for momentum diagonal 
fermion annihilation (creation). Furthermore, we use n 
for the particle density, 5 = 1 — n for the doping per 
site, z for the coordination number, and e\ — Ej^f/\A\ 
for the expectation value of the kinetic energy. The ket 
\M') = % fl-AO is the fully polarized Fermi sea of \- 
electrons from which one e^~ at the Fermi level ef is re- 



moved. The energy ba. 
infinite U {g = 0) (see * 



of 



with / = reads at 



Ae 



^/S-SF+EkSil-^/Szt) 2 ) 



(4) 



where Ek is the dispersion. The maximum energy low- 
ering is obviously obtained for k belonging to the lower 
band egde Eb, i.e. here fc b = 0. 

For finite / majority spin hopping processes from the 
position of the flipped spin to nearest neighbor sites are 
taken into account. This ansatz will be denoted NN. The 
amplitudes of these hopping processes are assumed to re- 
flect the lattice symmetry. Basile and Elser investigated 
an ansatz similar to NN which includes "f-hopping pro- 
cesses from the position of the |-electron to all other lat- 
tice sitestLy. Since the number of variational parameters 
increases with the lattice size they only studied a finite 
square lattice. The resolvent method developed in sect. [H] 
allows us to investigate a variational ansatz equivalent to 
the full Basile-Elser wave function in the thermodynamic 
limit on all common lattices. We also derive improved 
variational criteria for the Nagaoka instability at U < oo 
by extending the Hilbert subspace further. 



II. RESOLVENT APPROACH 

Generally, a resolvent is an operator-valued expression 
of the type a 



R{u>) = \/{u-{H-E N )) 



(5) 



where H — Ej^ is the Hamilton operator with respect to 
the ground state energy E^/ (here: the Nagaoka state 
energy). From (JsJ) it is clear that the existence of any 
state at lu = ljq implies a pole or at least a singularity 
in the resolvent. For this reason, we will investigate in 
the following the resolvent R applied to c\ , \Af') and 
compare u)q to e-p . 

It is not possible to compute R for the whole Hilbert 
space except under aLmplifying conditions like infinite co- 
ordination numbecj. Hence we will restrict the inversion 
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to certain subspaces which still allow an analytical treat- 
ment. The results obtained in this way for the lower band 
edge are variational. This means that excitation energies 
found are upper bounds to the true ones and that spe- 
cific interaction values U come out too small compared 
to those of the full solution. 



A. Case U = oo: Ansatz RESO 

For infinite on-site repulsion no double occupancy is al- 
lowed. Thus at the site of the |-e~ no f-e - is allowed. We 
investigate therefore the variational subspace spanned by 

a. 



^aJ^al^M') with arbitrary i and j. We define 



A-|AO 



Ak := | A| x ' 2 e Mi(k h - fc)l)%c+ ' 
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(6a) 



(6b) 



where the admissible values of k are outside the Fermi 
sphere (FS), but inside the Brillouin zone (BZ), i.e. 
k G BZ \ FS. The Hamiltonian does not mix states @ 
for different total momenta fc b . States (0) for different 
total momenta k h are orthogonal. Ansatz (JTJ) contains 
in particular the NN ansatz (0) and of course the sim- 
ple Gutzwiller ansatz. It comprises "["-hopping processes 
of arbitrary distance, i.e. it is the thermodynamic exten- 
sion pf. the ansatz investigated previously by Basile and 
ElseiH 

For the computation of the resolvent R(u>) one can 
use the Mori/Zwanzig projection formalism (see e.g. 
appendix C inEj) with the scalar product (A\B) := 
(Af'\[A + ,B] + \Af') for the operators A and B. The re- 
solvent d|) then becomes 



R 4l ^(w) = <*fc 1 |flM|$fc B > 

= (A hl \(u-C)- 1 A* a ) ■ 



(7) 



Here the Liouville operator L is ruscd which is defined as 
LA := [H, A] for all operatocsylES. The resolvent can be 
expressed in matrix notations by 

R(w) =P(cjP-L-M(cj))- 1 P (8) 

with the norm matrix P and the frequency matrix L 



\H-E N \<b K 



(9a) 



(9b) 



The frequency matrix L encodes the effect of H in the 
subspace considered. The deviation of P from unity ac- 
counts for the non-orthonormality of the basis. The so- 
called memory matrix M(u>) describes the effect of all 
processes which imply excursions outside the subspace 
considered. If the ground state is known exactly (which 



holds in the present case) the approximation M(oj) = 
is variational in nature for the lower band edge. 

It is the aim of the subsequent calculation to obtain 
a simple condition for the singularity of (uP — L). This 
singularity then signals that lu corresponds to an eigen 
energy. To this end, we first need the matrix elements 



= nS 



ii ±2 



|A|" 



Lfej ,fe 2 t 



fa.,kAn ■ £fc, - ei) 



£h - Sk i: k 2 (zty 



ei £ k 2 -k h 



(10a) 
(10b) 



(10c) 



We use the notation e, := (Q(sf — %)£|)bz (©(e) is the 
Heaviside function). The elements in (]!(]) are obtained 
with the help of Wick's theo rem si nce \AT') is a simple 
Slater determinant. In ( 10b ) and (|l0c|) , we distinguish 
the part coming from the motion of the "["-electrons and 
the part coming from the mot ion of the |-electron. The 
expression {zt)~~ x ei £fc 2 -fc in (10c) is obtained from 



-(2ir)- d J e k2 - hh - hi d d kx = 
fejeBZVFS 

Yt £ k 2 -k b ( 2 ^r d / e^dtki . (11) 



fc x GBZ\FS 



This relation holds for all homogeneous, isotropic lat- 
tices with NN hopping only, e.g. square lattice, triangu- 
lar lattice, kagome lattice and so on. The result (|ll| ) 
can be found easiest by interpreting the left hand side 
as convolution of £k -k and of 0(ef — Ek^) > i- e - as 
a multiplication in real space which concerns only the 
NN terms. Thus it is the multiplication with a constant 
(cij af) — —ei/(zt). The sites i and j are arbitrary ad- 
jacent sites since all bonds are equal due to the required 
homogeneity and spatial isotropy. 

On the basis of dlTj) the matrix inversion can be 
rephrased as 



(wP-L)- 1 = (d- 1 + (oJ-e h )vv+Y 



(12) 



with the constant vector v = |A| 1 / 2 and the diagonal 
matrix d*^ = 5/^/ fe) witn 

f(k) := [n(u - £fc) + ei(l + (zty'e^kjr 1 ■ (13) 

Note that the dyadic product vv + provides a S\A\ x 5|A| 
matrix with the constant matrix element |A| _1 . 

Expanding the right hand side of ( |l2"l) in terms of vv + 
and resummation in terms of 

h(oj) := v + dv_ 

= (27T)- d J f(k)d d k 
fc<EBZ\FS 

= (Q(e k - ef) /(A))bz (14a) 
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yields 

= (cT 1 + (w - e^vv+r 1 
= d(l + (uj-e h )vv + dr 1 

= d (1 — (w — £h)«f + d + (w — £h) 2 «« + dt;f + d + . . .) 
= d-(u-e b ) dvv + d (1 + (w - £ b ) • ( 14 b) 



The matrix elements thus read 
((wP-L)- 1 ), , =5 kid 



W - £b 



1 + (w - ffb)^(w) 



; 2 /fe) 
|A| 



(14c) 



From (14c) we read off that (cjP — L) is singular for 

= 1 + (w - £ b ) . (15) 

The trick to reduce dyadic perturbations to simple divi- 
sions is commonly known under the name 'Householder 
method' in the numerics of matrices. This extremely 
simple result is derived here for all Bravais lattices, e.g. 
the square lattice, the triangular lattice, but not for the 
honeycomb lattice or the kagome lattice. The restriction 
to Bravais lattices enters since we implicitly assume that 
there is one eigen state for each value of k in the one- 
particle Hamiltonian. But it will be shown in appendix 
|bJ that identical formulae apply for general unfrustrated 
lattices. Similar formulae can be found for frustrated 
non-Bravais lattices, for instance the kagome lattice in 
appendix |d| 

In appendix^ it is explained that computing h(u>) for 
non-bipartite lattices requires explicit integration over 
the momenta. For lattices where the band minimum £b 
is reached at k h = further simplification is possible. 
The band minimum is found at k h = in particular for 
bipartite lattices where one may choose t > without 
loss of generality. Then the term Ek-k b in ( ]l3| ) reduces 
t o th e unshifted dispersion and the whole integration in 
( |l4a| ) can be written as integration over the density of 
states (DOS) p(e). 



h(u>) 



P{£)d£ = i G(O) 



G(y) :- 



p(e) de 
y-e 



7 := n — ei/zt 
Q := (nuj + ei)/7 



(16a) 



(16b) 

(16c) 
(16d) 



We will call the ansatz deduced from the subspace 
give n in (|^) RESO. It leads to the singularity condition 
( |l5| ) or to its generalizations for non-Bravais lattices. 



Once the energy u> is found from (|15|) for a given Fermi 
energy £p the spin flip energy for the whole process of 
taking one |-e _ out at the Fermi level and inserting it as 
J,-e~ at the lowest possible energy is given by 



Aeoo = u - £ F 



(17) 



A critical doping <5 cr is found where this spin flip energy 
vanishes. 



B. Case U < oo: ansatzes RES1, RES2, and RES3 

Besides the calculation of variational upper bounds for 
spin flip energies and resulting critical dopings it is our 
aim to determine critical interaction values U. 

For U < oo we have to include states with double 
occupancy. The easiest way to do so is to include a local 
double occupancyoES. This is done in the ansatz RES1 
by adding to the states defined in <M) the state 



|*0 :=|A| 



-1/2 



^exp(ifc b i)a+o iT a+ 



m . (18) 



This ansatz contains the nearest neighbor ansatz NN (g) 
(and the Gutzwiller ansatz) for U < oo. Again we want 
to compute the resolvent (pi). To do so the parts com- 
puted in the previous subsection can be used again. The 
matrices for RES1 have the block structure 



Pi 


0+ 





P 2 


'wP 


L 



luP L = 



Di 


K 


N + 


D 2 



Bi 


M 




B 2 



(19a) 



(19b) 



The matrices Pi and Di arc the same as in (14c) at 
U = oo. The null vector in P comes from the fact 
that the state with double occupancy |*i) is orthogonal 
to the states without double occupancy |$fc). The other 
matrix elements are again found by Wick's theorem 

P 2 = n (20a) 

D 2 = n{u; -U) + e 1 -e b (n 2 - ( ei /(zt)) 2 ) (20b) 

N k = -lAp 1 / 2 (n(eb ~ £k) + ei(l + e k - k J{zt))) (20c) 

Since we are at present only interested in the singularity 
condition it is sufficient to compute one of the elements 
of oj~P — L. The easiest is B 2 , fo r which an argument 
similar to the one leading to ( |l4b| ), yields 

B 2 = (D 2 -N+B^ 1 N)- 1 . (21) 

Thus the singularity condition simply reads 

= D 2 - iV+Di 1 N (22) 
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Now it is advantageous that Di 1 is already given in 
(14c). Inserting ( |20| ) one obtains after some cancellations 



e b - nU(l + e h )h(uj)) = 



(23) 



Equation (^3|) is as simple as (JT^) and enables us to cal- 
culate critical U values explicitly. Setting uj = ef in (23), 
which according to ( [17| ) corresponds to vanishing spin flip 
energy, renders U cr directly accessible: 



U. 



RESl 



(S) 



£b 



(1 -8)[1 + (ep - £b) h(e F )] 



(24) 



It turns out, however, that the values for U CT from ( |24| ) 
are not very good close to half-filling n = 1 where anti- 
ferromagnetic exchange processes are important. These 
are not accounted for in fll8|). They are considered, at 
least to a certain extent, in the ansatz RES2 by using 



1*2 



:= IAI-V 3 £ expfeiK> jT a+|AO (25) 
<ii> 



instead of |*i) as extension of the RESO subspace. 

The block structure ( fL9[) re mains the same and so does 
the singularity condition (|22j) . Only the matrix elements 
are modified 



P 2 = (ei + 5e 2 )/t 2 

D 2 = {{e\ + 5e 2 ){uj -U) + e x e 2 + Se 3 

-£ b {e 2 1 ~e 1 e 3 /{zt) 2 ))/t 2 
Nk = 

|A|- 1/2 (ei(eb - £k) + e 2 (1 + £k-kj{zt))) /t 



(26a) 
(26b) 

(26c) 



The explicit expression resulting now from (|22| ) is less 
transparent than ( p3| ) since no cancellations occur. We 
focus here on the most important case k h — 0. In addi- 
tion to the definitions (Rq) we use 
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y 



ei - e 2 /{zt) 
(ei£ b + e 2 )/j' 

[ih) {8 + (n b -n)G(Q)) 



(27a) 
(27b) 
(27c) 



and obtain from (B2 



D2 = E ^( D i" 

i 



7 



-y 



1 + (lj - Eb)h(cj) 



(28) 



from which U cr can easil y be determined. The value U CI 
appears only in D 2 , see (|26b|) . The results of RES2 @ 
generically lead to J7 cr cx 1/5 on vanishing doping. In 
this sense it represents an important improvement over 
RESl d23]). For explicit results we refer the reader to the 
next section. 



At last in RES3, we generalize the variational states 
with double occupancy like ([l8]) and ( p5|) in the same 
manner as we generalized the states without double oc- 
cupancy in (||) 



B k \M') 



(29a) 



l A r V2 E ex P(*fe + k)i)a+c^a+ , (29b) 



where now the admissible values of k are all vectors inside 
the Fermi sphere (FS). Note that the extension RES3 
contains both RESl and RES2. The block structure of 
the resulting problem is similar to the one in (|l9|). The 
difference is that all blocks are now macroscopically large 



(30a) 



Pi 


0+ 





p 2 



(wP-L)" 



wP L 


■( 


D x 


N 




N+ 


D 2 


i ( Bx 


M 


■)■ 




- 1 V M+ 


B 2 







(30b) 



The matrix elements and details of the evaluation are 
given in the appendix The main problem is that one 
has to find a tractable condition for 



B 2 1 = D 2 N+Di : N 



(31) 



to be singular. But with expansion tricks similar to the 
ones used above this obstacle can be overcome. For bi- 
partite lattices a relatively simple final formula is found 
(|A22| ). An evaluation for the triangular lattice (appendix 
y) and the kagome lattice (appendix ^) is possible as 
well. 



III. RESULTS FOR VARIOUS LATTICES 



A. Square lattice 

The square lattice represents the simplest bipartite lat- 
tice structure in two space dimensions and has therefore 
been at the center of interest in most of the publica- 
tions dealing witti tbn variational investigation of Na- 
gaoka stability^^ElEjO'OEZlo. The energy band is given 
by 



£□ (fc) = — 2t (cos k x + cos k y ) , 



(32) 



with the lattice spacing set to 1. The DOS pn(e) which 
is depicted in fig. |l](a) can be expressed by a complete 
elliptic integral of the first kind (see appendix |]). For 
positive hopping matrix element t the lower band edge is 
reached at k h = while the maxima of the band structure 
are located at the corners of the square shaped first Bril- 
louin zone (fc t = (±7r,±7r)). The logarithmic van Hove 
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singularity at e = corresponds to the saddle points of 
the dispersion (^). The symmetric shape of the DOS 
with respect to e = reflects the particle-hole symmetry 
of the Hubbard model on the square lattice. In the fol- 
lowing we make use of this symmetry and consider only 
the case of a less than half filled lattice (0 < n < 1) and 
t > 0. 
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FIG. 1. (a) DOS for the square lattice (t = 1), (b) spin 
flip energy at U = oo as a function of the hole density for Gw, 
NN and RESO on the square lattice (t = 1). 



Fig. |l|(b) shows the spin flip energies at U — oo result- 
ing from the variational criteria discussed in the previous 
sections as a function of S. The Gutzwiller wave function 
((||) with / = 0) gives a critical hole density 5 CI = 0.4905 
for the instability of the Nagaoka stateEi For the varia- 
tional ansatz (|3|) including nearest neighbor hopping pro- 
cesses of the majority spins (finite /), the spin flip energy 
is considerably lowered and the critical hole density de- 
creases to 6 C1 - = 0.4155. The evaluation the variational 
state RESO, which contains all spin-up hopping terms 
of the Basile-Elser type, leads to <5 cr = 0.4045. Thereby 
we reproduce up to the fifth digit our result obtained 



inE-J where we took into account hopping processes over 
a distance of up to four lattice spacings. 

The fact that the reduction of the spin flip energy in fig. 
[j](b) is mainly due to the nearest neighbor term demon- 
strates the overwhelming importance of local polariza- 
tions of the spin up Fermi sea for the instability of the 
Nagaoka state. The resolvent method treats implicitly 
an infinite number of variational parameters and makes 
it 
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FIG. 2. (a) Phase diagram (n < 1): Nagaoka instabil- 
ity lines on the square lattice for Gw (dashed-dotted), NN 
(long dashed), RES1, RES2, RES3 (full lines, from hettom to 
top), and the 1100 parameter ansatz of Wurth et al.CJ (short 
dashed), (b) spin flip energy for U Ic d = 0.8 and t = 1 as 
a function of the hole density for Gw (dashed-dotted), NN 
(dashed), RES1, RES 2, and RES3 (full lines, from top to 
bottom). 

possible to investigate the full Basile-Elser ansatz for the 
first time in the thermodynamic limit. Compared to the 
iterative method used ino it has the remarkable advan- 
tage that the lowest possible spin flip energy in a given 
subspace can be calculated without explicit knowledge of 
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the corresponding state. As we will see in section E, it is 
not generally true that the best value for <5 cr within the 
Basile-Elser subspace can be obtained by restricting the 
spin up hopping processes to a small cluster centered at 
the position of the flipped spin. 

Fig. ||(a) shows the Nagaoka instability lines in the 
phase diagram for the Gutzwiller single spin flip (Gw), 
the nearest neighbor ansatz (NN) (^) as well as for the 
wave functions RES1, RES2, and RES3 evaluated by 
means of the resolvent method. The "("-hopping terms 
appear to be much less efficient in suppressing the Na- 
gaoka state if the hole density is small (because most of 
the sites near the flipped spin are already occupied by 
a spin-up electron) and the on-site repulsion U is finite 
(because the terms all exclude double occupancies at the 
down spin position). Since the Gutzwiller projector (with 
g > 0) represents the only term contained in RES1 which 
is relevant for U < oo, the critical on-site repulsion near 
half filling is only slightly increased and U cx remains fi- 
nite for 5 = 0. A remarkable improvement is obtained 
by allowing for nearest neighbor exchange processes and 
thereby taking into account the antiferromagnetic ten- 
dency of the nearly half-filled Hubbard model. This is 
embodied in the ansatz RES2. For a constant non-zero 
value of the DOS at the upper band edge it leads to the 
asymptotic behavior £/ cr ,red(<5) = 1—0(6) for S — > 0. This 
implies the instability of the Nagaoka state for all finite 
values of U in this limit. Fig. ||(b) shows that the opti- 
mum spin flip energy for RES2 plotted as a function of 
the hole density for a fixed finite value of U approaches 
a finite negative value of the order t 2 /U at half filling 
while it vanishes for all wave functions containing only 
the Gutzwiller projector. 

The asymptotic behavior for S —* of the spin flip 
energy and of the Nagaoka instability line U CI (8) is not 
affected by the extension of the Hilbert subspace to the 
full resolvent ansatz RES3. As for U = oo the local terms 
play the most important role in destabilizing Nagaoka 
ferromagnetism. With increasing hole density exchange 
processes become less important and the Nagaoka insta- 
bility lines for RES2 and RES3 approach the one ob- 
tained for RES1. Since all RES wave functions differ 
only in the subspace with double occupancies the corre- 
sponding instability lines end up with a diverging on-site 
repulsion U CI at the critical hole density 5 CT — 0.4045 
obtained for RESO. 

Fig. ^(a) displays also the best known variational 
bound for the Nagaoka stability regime on the square 
lattice computed by Wurth et alB3. The corresponding 
state contains 1100 terms, most of them describing exci- 
tations of the spin-up Fermi sea with up to two particle- 
hole pairs located within a 9x9 plaquette around the 
down spin position. The critical hole density obtained 
with this variational wave function is S cr = 0.2514 and 
the minimum critical on-site repulsion is U™ m /t = 77.74 
(RES3: U™ a /t = 36.21). Comparing these results one 
should keep in mind that the resolvent method allows 
to derive analytic expressions for the Nagaoka instability 



line U cr (S), at least for RES1 and RES2, while the cal- 
culation of the phase boundary for the 1100 parameter 
state requires an immense numerical effort. 



B. Square lattice with next-nearest neighbor 
hopping 

Extending the Hamiltonian ([[]) by taking next-nearest 
neighbor hopping processes of the electrons into account 
and introducing a corresponding hopping amplitude t' 
allows to create a particle-hole asymmetry of the DOS. 
Variation of the ratio t'/t makes it possible to simulate 
a continuous "transition" between a bipartite and a non- 
bipartite lattice. In this subsection we investigate how 
this transition affects the stability of the Nagaoka state 
with respect to a Gutzwiller single spin flip on the square 
lattice. Furthermore we will give a perturbation argu- 
ment for |i| < \t'\. 

The band dispersion of the so-called t-t'-U model on 
the square lattice is given by 

£t—v (t.) — —2t ( cos k x + cos k y ) — At' cos k x cos k y . (33) 



For t, t' > the lower band edge £b 



-A(t + t') is 



reached at k h = 0. The maxima of the band structure 
are located at the corners of the Brillouin square for t! < 
t/2 and at the edge centers for t' > t/2, respectively. 
Exactly for t' = t/2 the maximum single particle energy 
£t = 2t is reached at the whole border of the Brillouin 
zone. This leads to a nesting situation and to the largest 
possible partice-hole aymmetry with a diverging DOS at 
the upper band edge. For t' > t/2 local minima of the 
band structure develop at the corners of the Brillouin 
zone leading to a step in the DOS. In the limit t/t' — ► 
the single particle energy at these fc-points reaches the 
lower band edge. The calculation of the DOS p t -t'{e) 
requires in general a numerical fc-integration. Only for 
t' = t/2 it is possible to map pt-t' (e) on the DOS for t = 
and hence on a complete elliptic integral (see appendix 



ft-*(e) = (l-~ 



pu(2tjl-- 



(34) 



The symmetry of the Nagaoka stability regime with re- 
spect to half filling found in the "pure" Hubbard model 
is destroyed if the next-nearest neighbor hopping t is 
switched on. In analogy,-tp the non-bipartite triangular 
and kagome lattices (seec3 and sect. [II F in this paper) 
one should expect that the tendency towards saturated 
ferromagnetism increases for more than half filling and 
decreases for n < 1. The RES ansatzes with the reduc- 
tion to DOS integrals cannot be used for the t-t' model 
since the i'-hops go beyond nearest neighbor hopping. 

The calculation of the optimum spin flip energy for the 
Gutzwiller ansatz ((^) with / = 0) requires additional 
effort for the t-t'-U model due to the more complicated 
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structure of the band dispersion ( |33| ) . The kinetic energy 
of the nipped spin no longer depends only on eh but also 
on the corresponding momentum fc b . For t,i' > (i.e. 
for less than half filling) we find k h — as for if = 0, 
whereas for t, t' < (i.e. for more than half filling) we 
choose k h = (w,ir) for t'/t < 1/2 and k h = (tt,0) for 
t'/t > 1/2. 

Fig. I shows the DOS for the 



(a) 

1.0 n , 1 , r 




n 



FIG. 3. t-t'-U model on the square lattice for \t'\ < \t\/2 : 
\t'\ = 1 — |*| = 0, 0.1, 0.2, 0.3, 1/3: (a) DOS p t -t'(e) for 
t, t' > (pt_ t '(e) for t,t' < is obtained by e «-> — e), (b) 
Nagaoka instability lines for a Gutzwiller single spin flip (the 
curves for n < 1 correspond to t, if > whereas the curves 
for n > 1 correspond to t,t' < 0). 

t-t'-U model on the square lattice and the corresponding 
Nagaoka instability lines in the phase diagram for various 
ratios t'/t < 1/2. We set \t\ + \t'\ = 1 so that the lower 
band edge is always at £b = —4. Increasing t'/t leads to a 
lower DOS at £b and a higher DOS at e%, while the loga- 
rithmic singularity at e = At' approaches the upper band 
edge. The maximum particle hole asymmetry is reached 
at \t'/t\ = 1/2 (i.e. \t'\ = 1/3) where the DOS @ di- 



verges like (V £ t — £ • | log(e t — e )|) _1 for e ~ e t - The 
Nagaoka stability region for less than half filling shrinks 
as t'/t is increased and disappears at t'/t — 1/2 (Figs. 
||(b), ||). On the other hand, it expands rapidly for n > 1, 
especially in the limit n — 2. At t'/i = 1/2 the Nagaoka 
state is stable towards a Gutzwiller single spin flip for 
all t/ > in this limit. Even the slope of the Nagaoka 
instability line U CI (n) vanishes at n = 2. 

If one increases the ratio t'/t beyond 1/2, the loga- 
rithmic singularity in the DOS is gradually shifted back 
towards e = and the shape of pt-t' (s) becomes more 
and more symmetric (fig. ||(a)). 
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FIG. 4. t-t'-U model on the square lattice for \t'\ > \t\/2 : 
\t'\ = 1- \t\ = 1/3, 0.4, 0.5, 2/3, 1: (a) DOS p t -t'(e) for 
t,t' > (p t _ t '(e) for t, t' < is obtained by e <-» — g), (b) 
Nagaoka instability lines for a Gutzwiller single spin flip (the 
curves for n < 1 correspond to t, t' > whereas the curves 
for n > 1 correspond to t,t' < 0). 

Nevertheless the step at e = 4t'(l — t/t') remains present 
for all t/t' > 0. The DOS at t = is identical to p n (e), 
which reminds us that the t'-U model with suppressed 
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nearest neighbor hopping consists of two completely de- 
coupled square lattices. 

At t' = t the Nagaoka stability region in the phase di- 
agram is found to be still very asymmetric with respect 
to n = 1 (fig. |(b)). A further increase of t'/t makes the 
phase boundaries above and below half filling approach 
the ones obtained at t = 0. Within our variational calcu- 
lations, the local stability of the saturated ferromagnetic 
state is identical in both limiting cases t' = and t = 0, 
but see the perturbative argument below. The step in 
the DOS, however, leads to a cusp in the Nagaoka sta- 
bility line U cr (n) for all \t/t'\ < 2. For t — ► 0, this cusp 
approaches n = 1 and dU CI /dn\ n= i+ is discontinuous at 
t = 0. This represents a qualitative difference to the limit 

In fig. | the upper and lower critical densities for the 
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FIG. 5. Critical densities for the Nagaoka instability at 
U — oo for a Gutzwiller single spin flip on the square lattice as 
a function of \t'\ = 1 — Between the two full lines the Na- 
gaoka state is found to be possibly stable. The dashed-dotted 
line marks the singular case \t'/t\ = 1/2 where the parti- 
cle-hole asymmetry reaches its maximum. 

Nagaoka instability at U = oo are plotted as functions 
of \t'\ = 1 — \t\. Thereby we once again demonstrate 
the shift of the Nagaoka stability region towards more 
than half filling with increasing particle-hole asymmetry 
in the DOS. The regimes of complete Nagaoka stability 
for n > 1 (—0.21 < t' < —0.39) and of complete Nagaoka 
instability for n < 1 (1/3 < t' < 0.45) are not symmetric 
with respect to \t'\ = 1/3. There are two different reasons 
for this asymmetry. First, since the increase of the DOS 
at the lower band edge is more pronounced for t' \ —1/3 
(that is, on the left hand side of the dashed-dotted line 
in fig. g|) than for t' f —1/3, also the tendency towards 
saturated ferromagnetism in the low density limit (corre- 
sponding to n — > 2 in fig. |^) is stronger in the former case. 
Second, the Nagaoka instability condition near half fill- 
ing is essentially determined by the ratio e%/{zt), i.e. by 
the asymmetry of the band edges with respect to e = 0. 
The fact that the latter asymmetry is more pronounced 
for \t'\ > 1/3 than for \t'\ < 1/3 is responsible for the 



instability of the Nagaoka state for less than half filling 
on the right hand side of the dashed-dotted line in fig. |j. 

In the limit t — > 0, a perturbative arguments gives fur- 
ther insight in the stability of saturated ferromagnetism. 
Starting point is the observation that at t = the square 
lattice decomposes into two independent square lattices 
tilted by 45° with hopping element t'. Without any t 
the two independent Nagaoka states on each sub-lattice 
can be oriented arbitrarily without influencing the en- 
ergy. Thus we deal with a degenerate situation and in- 
vestigate by (second order perturbation coefficient 
in t) whether the parallel or the antiparallel orientation is 
favored. The linear order vanishes for particle-hole 
symmetry reasons and does not lift the degeneracy. 

For the parallel configuration it is straightforward to 
calculate E^ 2 \ Without loss of generality we choose t' = 
1/4 and consider e{k) = Eo(k) — 2t(cos(k x ) ± cos(k y )) 
with £o(k) = — cos(fc x ) cos(fc ?/ ) as dispersion. The plus 
sign refers to n < 1 , t' > and the minus sign to n > 
l,t' > 0. This can be seen by means of a particle-hole 
transformation and a sign transformation Cj — * — a on 
all sites with even x-coordinate. One obtains at constant 
filling = - \A\/(2t')A±(e F ) with 

A±(e F ) = 

/d 2 k 
(cos(fc x ) ± cos(fc y )) 2 (5(e F + cos(k x ) cos(k y )) 

— TV 

= —$ s (±e F K(l-4)-E(l-4)) (35) 
yielding the dotted curves in fig. ^. The relation 

i 1 1 r ,- 1 1 1 , 1— 
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n 
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1 .4 



FIG. 6. Second order perturbation coefficient in t 

in units of At' N . Dotted line: for parallel Nagaoka states 
(or global ferromagnetic state, see text); Dashed line: up- 
per bound to E' 2 ' for antiparallel Nagaoka states (or global 
antiferromagnetic state). 

( p5| ) is found with the help of the quantities /„ in ap- 
pendix A of Hanisch/Mullcr-HartmanrO; K and E are 
complete elliptic integrals. Note that the coefficient E^ 
is not continuous across n = 1. 
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Next we assess the energy of two antiparallel Nagaoka square lattice at t 



states on each of the sub-lattices. Let us use a\ a for the 
fermions on the A sub-lattice and b\ a for the fermions 
on the B sub-lattice. The perturbation reads then 

Hi = -2t ( C0S (M ± C0S ( k y))( a t* b k,v + b t,a a k.a) 



as no saturated ferromagnetic 



fceMBZ 



(36) 



where MBZ is the magnetic Brillouin zone. The second 
order energy lowering is 

E {2 h 2 \k\ = -(A],B I ^{Ho-Eq^H^A^B I)) . 

(37) 

The acronyms A | and B J, stand for the respective Fermi 
seas. There are two processes which contribute equally 
to (|37]) . Either a fermion is shifted from A to B and back 
or a fermion is shifted from B to A and back. The latter 
yields explicitly 

£(2)=_ m 51 (cos(M±cos(fc. y )) 2 e(e F -£oCfc) 

' ' fcGMBZ 

(A t \a k>1 {H 0tA - E ,a ~ e©)"^^! A t) (38) 
d 2 k 

———(cos(k x ) ± cos(k y )) 2 g k (e (k)Q(e F ~ e (k) 

where gu is the one-particle Green function. Now we 
specify that we work at U = oo and we assume that the 
Nagaoka state is stable for t — at the filling considered. 
If the Nagaoka state is not stable we do not need to make 
the present comparison anyway. Based on our assump- 
tion, the Green function is purely real and negative. It 
obeys the inequality 



SfcMfc)) < MAO-efe)- 1 < 

efc := (a^ ^Ho^A - E 0i a\oI^) 
= -e l /5 + e a {k)5{\~(ei/5f) 



(39a) 



(39b) 



The estimate ( |39a| ) corre sponds to a simple Gutzwiller 
ansatzES and yields ( |39b| ) (see (4) and (5) with t = t' = 
1 /4 and z = 4 inO) . Thus we obtain 



E m < « 



ds 



A ± (e) 



(40) 



where / = S(5-l)-ef, A = ex/f, and A± from @. The 
evaluation of the right hand side of (M) yields the dashed 
curves in fig. ||. The essence of fig. |gis that the saturated 
ferromagnetic state is unstable in the limit t — ► for all 

(2) 

fillings. The small region where Ep-^ lies below the upper 

(2) 

bound for E A p M does not count since we know that at 
these dopings (and for larger dopings) already the pure 



ground state, see e.g.E3B 

We wish to draw the reader's attention to the fact that 
the comparison in fig. || is quite different from the main 
theme of this paper which is based on single spin flip 
energies. Here the global stability is tested with a com- 
pletely different, antiferromagnetic state. We learn from 
the perturbative argument that in fig. [| the true lines 
n{t') comprising the global Nagaoka stability region have 
to converge both to the point t' = l,n = 1. 



C. Simple cubic lattice 



The energy dispersion of the simple cubic lattice is 



e sc (k) = — 2t (cos k x + cos k y + cos k z ) 



(41) 



The calculation of p sc (e) can be performed by an inte- 
gration over the known DOS of the square lattice (see 
appendix |Ej). The maxima and minima of the energy 
dispersion (|4l[) are Et = z\t\ and Sh — — z\t\, respectively, 
with the coordination number z = 6. At the band edges 
the DOS (fig. 0) shows the square root 



. 15 




FIG. 7. DOS for the simple cubic lattice (t = 1). 

behavior which is characteristic for d = 3. The van Hove 
singularities at e — ±2t correspond to the saddle points 
of the dispersion ( |4l| ) . 

Fig. ||(a) shows the spin flip energy at U — oo for Gw, 
NN, and RESO. For small hole doping the loss of spin-up 
kinetic energy due to the spin flip is sufficiently strong 
to keep the Nagaoka state stable. With increasing (5 the 
spin flip energy decreases due to the gain of kinetic energy 
for the flipped spin which grows linear with 5 in leading 
order. The upper bound for the, critical hole density is 
reduced from S cr = 0.323 for GvM<B to 5 CI = 0.247 for 
NN and finally to <5 cr = 0.237 for RESO. As in d = 2, the 
NN hopping term gives the dominant contribution to the 
decrease of S cr while the extension of the spin-up hopping 
processes to the whole lattice has only a small effect. 
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FIG. 8. (a) Spin flip energy at U — oo as a function of the 
hole density on the sc lattice (t = 1) for Gw, NN, and RESO, 
(b) phase diagram (n < 1): Nagaoka instability lines on the 
sc lattice for Gw (dashed-dotted) , NN (long dashed), RES1 
(lower full line), RES 2 (short dashed), and RES3 (upper full 
line) . 

RothH investigated the Nagaoka instability with re- 
spect to a single spin flip on the sc lattice already in 1969, 
making use of the so-called two pole approximation ip-. 
stead of the projection method. It was shpwn latero'tHl 
that the Hilbert subspace considered ir£3 is equivalent 
to the Basilc-Elser subspace in the limit U — > oo. Roth 
obtained numerically a critical hole density of 0.24 which 
is consistent with our variational result for RESO. 

The phase diagram (fig. ||(b)) for the simple cubic lat- 
tice shows a qualitative difference to the square lattice: 
The critical U at half filling obtained for the Gutzwiller 
single spin flip is not at all improved by including NN 
hopping terms. Even for RES1 U cr (6 = 0) is still given 
by the band width 12|t|. This is due to the fact that for 
the sc lattice the DOS at the upper band edge vanishes 
while it is nonzero for the square lattice. 

As in d = 2, the ansatz RES2 leads to U cr (S = 0) = oo 



and to a considerable reduction of the Nagaoka stabil- 
ity regime near half filling. For the full resolvent ansatz 
RES3 we finally achieve a minimum critical coupling of 
U™ in = 48.9|i| (corresponding to U rcd = 0.753) below 
which the Nagaoka state is proven to be unstable for all 
S. 
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FIG. 9. (a) DOS for the bec lattice (t = 1). (b) phase 
diagram (n < 1): Nagaoka instability lines on the bec lattice 
for RES1 (dotted line), RES2 (dashed), and RES3 (full line). 

The region left for a possible Nagaoka ground state 
on the sc lattice is therefore substantially smaller than 
on the square lattice (RES3 for the square lattice: 
S cr = 0.405, C/ C n r lin = 36.2|i|). Generally the tendency 
of the Hubbard model towards a saturated ferromag- 
netic ground state on a d-dimensional hypercubic lattice 
becomes weaker with increasing d. Muller-HartmannLj 
showed that the critical hole density at U — oo with re- 
spect to a Gutzwiller single spin flip decreases asymptot- 
ically as S CI oc 1/yd In d for d ^> 1. In the limiting case 
of infinite dimensions the ground state of the Hubbard 
model is never fully polarizedcll. 
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D. bcc lattice 

The bcc lattice is another interesting example of a 
three-dimensional bipartite lattice. It has a slightly 
higher coordination number z — 8 compared to the sim- 
ple cubic lattice. Its dispersion reads 



£hon (k) = ± J t (3* - E A (fc) ) , 



(43) 



^bcc(k) = — 8tcos(fe x ) cos(ky) cos(fc 2 



(42) 



The calculation of the DOS Pf, cc (e) can again be per- 
formed by an integration over the known DOS of the 
square lattice (see appendix [^) . The bipartiteness is ob- 
vious since ebcc{k + Q) + Sbcc(k) — with Q = (it, ir, it). 
For this reason we consider only n < 1. 

The DOS is shown in fig. |^(a). The square root sin- 
gularities at the band edges are generic for three dimen- 
sions. The least common feature for a three dimensional 
lattice is the squared logarithmic singularity at zero en- 
ergy pbcc(s) ~ ln 2 (e)/(47r 3 ) which results from the points 
in momentum space where all cosines in ( (42| ) vanish, e.g. 
£bcc(k) « -8t(k x - n/2)(k y - %/2){k z - tt/2). 

Evaluating ([I]) for the bcc lattice, we find the critical 
density <5 cr = 0.324 in the Gutzwiller approach. This is 
almost the same result as for the simple cubic lattice. 
The result S cr = 0.239 for the full ansatz RESO is also 
only a tiny bit higher than the RESO critical doping for 
the sc lattice. It appears that the essential ingredients 
are indeed the dimensionality and the bipartiteness as we 
will see below. 

The results for finite interaction are shown in fig. ^(b). 
The value of C/br is 16.413. The reduced interaction 
values are very similar to the ones for the simple cubic 
lattice. The ansatz RES1 does not capture the diverg- 
ing interaction for n — > 1 but RES2 yields already the 
asymptotic behavior of RES3 for n — > 1. The critical 
interaction is C/ rc d, C r = 0.7438 for RES3. 

As far as the local stability of the Nagaoka state is 
concerned we do not find any indication that the bcc 
lattice is more favorably lhan the simple cubic lattice. 
Herrmann and NoltingEjO found in the framework of 
the spectral density approach an enhanced tendency to- 
wards ferromagnetism for the bcc lattice. They inves- 
tigated the divergence of the susceptibility in the para- 
magnetic phase which is enhanced by the large DOS at 
zero energy. Combining their result with ours one might 
come to the conclusion that the bcc lattice favors a non 
saturated ferromagnetism for intermediate coupling and 
doping. 



E. Honeycomb lattice 

Besides, the square lattice the honeycomb lattice (see 
fig. 9 ir£j) is another prominent example of a bipartite 
lattice in d = 2. In contrast to the square lattice it is not 
a Bravais lattice, however, but a triangular lattice with 
a two site basis. The coordination number is z = 3 and 
the band dispersion reads 



where £a(^) stands for the energy dispersion of the tri- 
angular lattice to be described in (|44|). Despite this ad- 
ditional complication the formulae developed in sect. [□] 
via the resolvent method hold here as well (see appendix 
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FIG. 10. Spin flip energy at U = oo as a function of the 
hole density on the honeycomb lattice (t = 1) for Gw, NN, 
and RESO. 



e with respect to Gw 
. Here we present the 



The instability of the Nagaoka st 
and NN was already discussed int 
improvements obtained by the resolvent method. The 
evaluation of RESO shows that hopping processes with a 
larger distance from the down spin position have only a 
very small influence on the optimum spin flip energy at 
U = oo (fig. £3). 
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FIG. 11. Phase diagram (n < 1): Nagaoka instability lines 
on the honeycomb lattice for Gw (dashed-dotted), NN (long 
dashed, almost identical with RES1), RES2 (short dashed), 
and RES3 (full line). 

The instability gap (0.379 < 5 < 0.481) between the 
two possible Nagaoka stability regions remains almost 
unchanged compared to the result for NN. The upper 
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critical hole density is only slightly improved to S CI = 
0.643 from 0.662 (NN) and 0.802 (Gw)El 

As explained inEa the Nagaoka stability island in the 
phase diagram around quarter filling (fig. is mainly 
due to the zero in the DOS at e = 0, i.e. between the 
two energy bands. Since the lattice structure enters the 
calculation of the optimum spin flip energy by means of 
the resolvent method only via the DOS the stability is- 
land is present even for the full resolvent ansatz RES3. 
On the other hand, the critical U at half filling diverges 
for RES2 and RES3 and the Nagaoka stability region for 
small S shrinks compared to the results for NN and RES1. 
These results and the pronounced difference between the 
two minimum values of U (41.2 \t\ for the low doping 
regime and 17.25 \t\ forr-the stability island) corroborate 
the previous conjecturec3 that a saturated ferromagnetic 
ground state exists around quarter filling. The Jack of 
a Nagaoka theorem for the honeycomb latticetjO indi- 
cates a degeneracy between the Nagaoka state and other 
possible states near half filling even at U = oo. 

F. Triangular lattice 



proving the instability of the Nagaoka state in the low 
density limit. The resolvent ansatz RES0 lowers the spin 
flip energy further and implies 5 CT = 0.824 (fig. |l2|(a)). 
The difference Ad = 0.088 between the results obtained 
for NN and RES0 is eight times larger than the one for 
the square lattice (A<5 = 0.011). 
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The triangular lattice is non-bipartite. It can be de- 
composed into three sub-lattices, each of them having 
triangular structure. Investigating the local instability of 
the Nagaoka state towards a Gutzwiller single spin flip 
a Nagaoka ground state was excluded on the triangular 
lattice for less than half fillingocj. This is in agreement 
with the Nagaoka theorem, which predicts a saturated 
ferromagnetic ground state at U = oo only for the half 
filled lattice plus an additional electron. Thus we con- 
sider henceforth the electron doped case for t = 1 or, 
equivalently, t = — 1 and n < 1. 

Each lattice site has z = 6 nearest neighbors located 
at the corners of a hexagon. The band dispersion is given 
by 

£A (fc) = 2cosfe x +4co S ^cosf^^ (44) 

where k belongs to the likewise hexagon-shaped first Bril- 
louin zone. The upper band edge (et = 6) is found at 
the center of the Brillouin zone, whereas the lower band 
edge £b = — 3 is reached at the corners_pf the hexagon. 
The DOS (see appendix |] and fig. 1 ir£3), which can be 
expressed by a complete elliptic integral, displays a log- 
arithmic van Hove singularity at e = —2. For ep = — 2 
the Fermi surface forms a hexagon with an area of 3/4 of 
the whole Brillouin zone. As usual in d = 2 the DOS at 
the band edges is nonzero (pb = Apt = (VStt)^ 1 ). 

In contrast to the square lattice, the Nagaoka state 
repahis stable towards Gw for all fillings n > 1 at U — 
ooc3i23. The corresponding spin flip energy as a function 
of S is depicted in fig. |l^(a). Evaluating NN, however, 
a negative spin flip energy is found above S cr = 0.912 



(b) 




6 = n - 1 



FIG. 12. (a) Spin flip energy at U = oo as a function of the 
hole density on the triangular lattice (t = — 1) for Gw, NN, 
and RES0, (b) spin flip energy at f7 re d = 0.4 as a function of 
the hole density on the triangular lattice (t — —1) for RES1 
(dashed-dotted line), RES2 (upper full line), and RES3 (lower 
full line). 

This demonstrates the importance of the spin- up hopping 
processes for the instability of Nagaoka ferromagnetism 
on the triangular lattice. The reason is that due to the 
large hole densities under consideration, the probability 
to find unoccupied sites near the flipped spin is quite 
high. The same line of reasoning applies also for U < oo, 
see fig. |l3[ 
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FIG. 13. Phase diagram (n > 1): Nagaoka instability 
lines on the triangular lattice for Gw (dashed-dotted) , NN 
(lower full line), RES1 (short dashed), RES 2 (long dashed), 
and RES3 (upper full line). 

Previously we investigated a variational state which re- 
stricts the hopping processes to a_31-site cluster around 
the position of the nipped spinO and obtained 5 CI — 
0.887. Although the nearest neighbor processes once 
again are the most important ones, the number of rel- 
evant hopping processes on the triangular lattice turns 
out to be much larger than on the square lattice. Hence 
the critical hole density 5 CI — 0.824 found by the resol- 
vent method is essentially lower than the one found from 
the finite cluster calculations. Moreover, the evaluation 
of RESO requires much less analytical and numerical ef- 
fort than the iterative extension of the variational ansatz 
by additional hopping processes. For details on the ap- 
plication of the resolvent method to the triangular lattice 
see appendix |C]. 

Near half filling the influence of the majority spin hop- 
ping processes contained in NN and RES1 (which sup- 
press double occupancies) on the Nagaoka stability is 
negligible, as expected (fig. [l3|). In contrast to this the re- 
solvent ansatz RES2 with nearest neighbor hopping pro- 
cesses creating double occupancies leads to a negative 
spin flip energy near half filling for all C/ < oo and hence 
to a divergence of U CT (n = 1) (fig. |l3|). It turns out, how- 
ever, that for larger hole densities, when the exchange 
effect looses its importance, RES2 is somewhat less suc- 
cessful than RES1. The plot of the spin flip energy as 
a function of S for the comparatively small on-site re- 
pulsion [/red = 0.4 in fig. [l2](b) demonstrates that above 
5 ~ 0.12 the creation of extra holes near the flipped spin 
as described by RES2 is energetically unfavorable. The 
full resolvent ansatz RES3, comprising RES1 and RES2, 
gives of course the best lower bound for the Nagaoka 
instability line U cr (6). The minimum critical coupling 
obtained for RES3 is [/™ in = 9.62|t| (C7 red = 0.378), the 
critical hole density at U = oo is given by the RESO 
value 6 CT = 0.824. Hence the region for a possible Na- 
gaoka ground state on the triangular lattice appears to 
be much larger than on the bipartite square and honey- 



comb lattices. 



G. Kagome lattice 

Taking the kagome lattice as an example of a frus- 
trated non-Bravais lattice we want to demonstrate that 
the resolvent method works .also for this class of lattices. 
Representing the line grapho of the honeycomb lattice 
the kagome lattice (for t < 0) shows a flat, i. e. disper- 
sionless band with spectral weight 1 /3 at the lower band 
edgee b = -2|t| (fig.0). 




FIG. 14. DOS for the kagome lattice (t = -1). 

All line graphs display such a flat bandEI. The kagome 
lattice is the first and the most neopiinent example of 
so-called flat- band ferromagnetisraiirca. A macroscopic 
degeneracy of the lowest single particle energy leads for 
certain band fillings to. a unique saturated ferromagnetic 
ground state. MielkeCJ proved that the Nagaoka state is 
the unique ground state of the Hubbard model on the 
kagome lattice for all U > at n = 1/3. Although in the 
flat-band regime every ground state of the Hamiltonian 
(0) is a simultaneous eigenstate of H^- m and H pot , the 
uniqueness of the ground state is not trivial. Eor n < 1/3 
the fully polarized ground state is not uniqueo. 

The kagome lattice can be considered asr-a triangu- 
lar lattice with a basis of three lattice pointsE3, see also 
appendix Besides the flat band e(k) = 2t the diago- 
nalization of -ffkm leads to the two dispersive bands 



e K (k) = -t[l±J3-£ A (k)/t 



(45) 



where e a {k) stands for the dispersion ( f44| ) of the triangu- 
lar lattice. For the kagome lattice the resolvent method 
requires less effort than for the triangular lattice with 
t < since the lower band edge £b = — 2|£| is reached 
at fc b = for one dispersive band and for the flat band 
of course. Thus q = is the optimum momentum as for 
bipartite lattices. Hence all lattice dependent quantities 
appearing in our formulae can be calculated as integrals 
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over the DOS p K {s)- For t = -1 one finds the DOS of 
the kagome lattice (see fig. |l^ and appendix |e|) as 

PK (e) = i 6(e + 2) + | e - 1| • M ((e - l) 2 - 3) ; (46) 

Pa(e) is the DOS of the triangular lattice. Nevertheless, 
the fact that the kagome lattice is not a Bravais lattice 
induces some changes in the analytic expressions for the 
spin flip energy (appendix |d|) . 
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FIG. 15. Spin flip energy at U = oo as a function of the 
hole density on the kagome lattice for Gw (dashed-dotted) , 
NN (dashed), and RESO (a) for t = 1, (b) for t = -1. 

Fig. |l5|(b) shows the spin flip energy for U = oo and 
t = —1 as a function of 5 foc-JlESO compared to the 
Gw and NN results obtained inc3. As for the honeycomb 
lattice the enhancement of the Nagaoka stability for 5 — > 
1/3 is due to the zero in the DOS. The effect of the 
additional spin-up hopping processes contained in RESO 
is most pronounced for S > 1/3. But the spin flip energy 
remains positive for all band fillings. Note that the exact 
result in the flat-band regime is a zero spin flip energyuO. 
The phase diagram for n > 1 in fig. [l6| 
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FIG. 16. Phase diagram: Nagaoka instability lines on the 
kagome lattice for Gw, NN, RES1 (lower full line), RES2, 
and RES3 (upper full line). For n < 1 the difference between 
RES1, RES2, and RES3 is less than the line width. 

shows a strong tendency towards Nagaoka ferromag- 
netism also beyond the flat-band regime, where we find 
the Nagaoka state to be stable for all U > 0. There is 
only a marginal difference between the Nagaoka insta- 
bility lines for NN and for RES1, since the values of U 
under consideration are too small to allow a significant 
reduction of the spin flip energy by Basile-Elser hopping 
processes. Near half filling, however, we are able to re- 
strict the Nagaoka stability region by RES2, i.e. by tak- 
ing antifcrromagnetic exchange processes into account. 
As for the triangular lattice, for a certain range of fill- 
ing around n = 3/2 away from half filling RES2 gives a 
weaker bound for U CI (8) than RES1. 

For positive hopping matrix element t the flat band is 
found at the upper band edge. The flat-band regime for 
n < 1 corresponds to hole densities < 6 < 1/3. Since 
fc b = and Sb = —zt the resolvent method formulae 
are those of the bipartite lattices (see appendix^). Figs. 
|l5|(a) and |l6| show-that the small Nagaoka stability island 
found previouslyo for very large U around quarter filling 
is still present for RESO - RES3. The upper critical hole 
density is reduced from 0.727 (NN) to 0.715 (RESO) and 
C/ C n r lin reaches 191.5|t| (RES3) instead of 129.5|t| (NN) 
though. These results may indicate that this stability is- 
land really provides an example of a saturated ferromag- 
netic ground state on-a non-bipartite lattice for less than 
half filling. Its origincS is the zero in the DOS enhancing 
the stability of the Nagaoka state around S = 2/3. 

H. fee and hep lattices 

The fee and hep lattices as the most prominent close- 
packed lattices in d = 3 are found in numerous real sub- 
stances among them the ferromagnetic transition metals 
Ni (fee) and Co (hep). The face-centered cubic lattice is 
a Bravais lattice with coordination number z = 12. Its 
band dispersion 
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— At (cos k x cos k y 



is related to the dispersion 
lattice viao 



cos k x cos A: z + cos k y cos fe z ) (47) 
) of the simple cubic (sc) 



Sfcc(k) 



It 



:(2fc) 



(48) 



The hexagonal close-packed lattice (also with z = 12) is 
not a Bravais lattice but a hexagonal lattice with basis. 
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FIG. 17. (a) Identical DOS for the fee and hep lat- 
tices (i = —1) and DOS for the fee lattice with additional 
next nearest neighbor hopping (t' — t/2), (b) phase diagram 
(n > 1): Nagaoka instability lines on the fee lattice for Gw 
and NN and on the t-t'-fee lattice with t' = t/2 for Gw. 

Within the hexagonal planes, which we assume to be par- 
allel to the xy plane, the fc-dependence of the two energy 
bands reduces to the energy dispersion (44) of the trian- 
gular latticeEj: 

£hcp(]£) — SA^X; ky) 

± 2tcos (V273fc 2 ) • ^3- e A {k x ,k y )/t . (49) 



Hence it is possible to compute the DOS of the hep lattice 
by integration over p&(e) (see appendix ||). 

Modelling the fee and hep structures with close-packed 
spheres, the sequence of layers with different positions of 
the sphere centers is known to be ABABAB... for the 
hep and ABCABC... for the fee lattice. Diagonalizing 
the kinetic part of the Hamiltonian ([[]) in each of the 
hexagonal planes it turns out that the terms reflecting 
the different arrangement of the planes disappear if one 
chooses the Fourier transformation in a convenient way. 
The densities of states for the fee and for the hep lattices 
are therefore identical as well as our variational results 
on the stability of the Nagaoka state with respect to Gw 
and NN. 

For less than half filling a saturated ferromagnetic 
ground state was already excluded due to the complete 
instability of the Nagaoka state towards Gw at U — oota. 
Therefore we only investigate the case of more than half 
filling which corresponds to t < 0. The DOS (fig. 0(a)) 
displays the d — 3 square root behavior p(e) oc y^t — £ 
at the upper band edge e t = z\t\. For the fee lattice 
the lower band edge = — 4|i| is reached on different 
lines in fc-space which intersect in several critical points 
located on the border of the Brillouin zone. This reduces 
the "effective dimensionality" of the van Hove singular- 
ity by one and leads to a logarithmic singularity in the 
DOS. As a consequence, the Nagaoka state remains sta- 
ble in the low density limit (corresponding to n — * 2 
for t > 0) for all U > with respect to Gw and NN (fig. 
|l7|(b)). This result indicates the strong tendency towards 
Nagaoka ferromagnetism on the fee lattice, especially in 
comparison with the triangular lattice where we proved 
the instability of the Nagaoka state even for U = oo in 
the low density limit. Also for low and intermediate hole 
doping the extension of the Gutzwiller wave function by 
nearest neighbor hopping processes yields only a slight 
reduction of the Nagaoka stability region in the phase 
diagram. This is in sharp contrast, for example, to the 
situation on the sc lattice. The resolvent method was 
not applied to the fee and hep lattices since three dimen- 
sional momentum integrals would have to be performed 
in order to calculate h(uj) and h(uj) (s ee ap pendix Q). 

As for the square lattice (see section [HI B ) the particle- 
hole asymmetry of the DOS is even enhanced if one 
extends H^ in by electron hopping between next nearest 
neighbor sites with a hopping amplitude t' . On the fee 
lattice, these sites form a simple cubic structure such that 
the additional contribution to the dispersion exactly com- 
pensates the second term on the right hand side of ( |4^ ) 
if t' = t/2. In this case the DOS for the t-t'-fee lattice is 
connected to the DOS of the sc lattice via 



3(1 -i) 



Psc 



(50) 



and therefore finally simplifies to an integral over p n (e) 
(see appendix ||]). The next nearest neighbor hopping 
with amplitude t' = t/2 creates a square root divergence 



1G 



of the DOS at the lower band edge (for t, t' < 0) in 
contrast to the logarithmic singularity obtained for t' = 
(see fig. 0(a)). ' 

The Nagaoka instability line for a Gutzwiller single 
spin flip on the fee lattice with t' = t/2 is compared 
with the result for the simple fee lattice (t! = 0) in fig. 
0(b). The more pronounced singularity of the DOS at 
the lower band edge leads to an even more pronounced 
stability of the Nagaoka state in the low density limit. 
We find U CI (S) ocl-S instead of U CI (S) oc l/log(l - 6) 
for t' = 0. The slight increase of the critical U at half 
filling is due to the different band width of the t-t'-U 
model (18|i| instead of 16|i|). 



IV. CONCLUSIONS 

In summary, we investigated the stability of the Na- 
gaoka state for a series of two- and three-dimensional 
lattices: the square t-, the square t-t'-, the simple cubic, 
the bee, the honeycomb, the triangular, the kagome, and 
the fee (hep) lattice. The results were mostly variational 
in nature and concerned the energy change due to a single 
spin flip. By the resolvent approach the eigenvalue prob- 
lem in the variational subspace was reduced to a matrix 
inversion problem. The relatively simple structure of the 
matrices under consideration permits to convert the ma- 
trix inversion into a scalar inversion (or the inversion of 
a 2 x 2 matrix). For the t-t' square lattice a perturba- 
tive approach in t was used as well for investigating the 
t/f -> limit. 

The ansatzes RESO - RES3 are particularly simple 
for unfrustrated, isotropic, homogeneous lattices with 
nearest neighbor hopping (sect. |J, appendix ^). For 
frustrated, non-bipartite lattices our approach is still 
tractable, though more cumbersome. To demonstrate its 
tractability we derived formulae for the triangular lattice 
(non-bipartite, Bravais lattice) and for the kagome lattice 
(non-bipartite, non-Bravais lattice). 

We believe that our variational criteria are well suited 
in order to investigate the lattice dependence of satu- 
rated ferromagnetism in the Hubbard model since they 
cover the most relevant local excitations in the spin-up 
Fermi sea but are still simple enough to be evaluated 
routinely on various lattices in various dimensions. Lo- 
cal and non-local band narrowing effects are present in 
these approaches. Since this fact is not obvious in the 
complete approaches we resort to the previous result (Q) . 
The factor 6(1 — (ei/Szt) 2 ) clearly describes the band 
narrowing of the flipped spin. It comprises two factors 
one of which is local (S) and thus survives also in the 
limit d — ► oo. The other factor (1 — (ei/6zt) 2 ) is very im- 
portant as well since it vanishes equally_pn 8 — > 0. But 
in the limit z — > oo on scaling t oc l/y^SiS the latter fac- 
tor degenerates to unity. This clearly shows its non-local 
character. Note the importance of the sequence of limits. 
The more sophisticated variational approaches discussed 



in the present work (RESO - RES3) comprise the ansatz 
(||). Thus they contain also local and non-local band 
narrowing effects. The other main effect is a direct en- 
ergy lifting of the minority electron due to the infinite 
(or large) on-site repulsion. Since the minority electron 
blocks a site, the majority electrons loose the kinetic en- 
ergy related to hopping onto or from this site, namel y e\ . 
This is seen best in the kinetic matrix elements in ( |10b| ) 
or in the energy denominator (0). Futhermore, we like 
to draw the reader's attention to the non-orthogonality 
as it can be discerned in (10a). It is very difficult to 



comprehend its effect intuitively. But we know from the 
extensive efforts to reduce the critical doping by including 
more and more correlations!^! that this non-orthogonality 
hinders the spin flip to gain enough energy to destabilize 
the Nagaoka state. The added states do not reduce the 
critical doping any further since they do not really en- 
hance the accessible Hilbert space. 

Besides the achievement of easily evaluated ansatzes 
the comparison of the phase diagrams presented here 
yields the following main results. For bipartite lattices 
the possible Nagaoka region shrinks rapidly with increas- 
ing coordination number z (cf. square and simple cubic 
lattice). Herrmann and Nolting did not investigate low- 
dimensional lattices because they suppose that ferromag- 
netism is-ejtcluded in d = 1 and 2 by the Mermin- Wagner 
theoremtaO. Note, however, that neither the Mermiit. 
Wagner theorem makes any statement on ground statescll 
nor any extended theorem can exclude a ferromagnetic 
ground state since the total spin as conserved quantity is 
not affected by quantum fluctuations. The shrinking of 
the Nagaoka region on increasing coordination number 
can be understood from the lowering of the DOS at the 
band edges or, equivalently, as effect of a lower and lower 
band edge £b- 

For the above reasons we investigated low-dimensional 
non-bipartite lattices where low DOS at the lower band 
edge can be avoided. Indeed, we found that the possible 
Nagaoka regions are enlarged considerably. This is true 
for electron doping for conventional hopping (t > 0, n > 
1) whereas saturated ferromagnetism in the hole doped 
region (t > 0, n < 1) can be excluded by our results. 
Treating the electron doping also as hole doping after 
a particle-hole transformation, i.e. t > 0, n > 1 — > t < 
0,n < 1, this phenomenon is easily understood: e^it < 
0) < Eb(t > 0). The ratio of the lower band edges is 2 for 
the triangular and the kagome lattice, and 3 for the fee 
and hep lattice. In infinite dimensions it b ecomes even 
oo for the generalizations of the fee latticeE3'L§EJ. For 
these generalizations one has £b(t > 0)/eb(t < 0) oc \fd. 

The above observations concern already the asymme- 
try of the density of states. Our results clearly show that 
a large asymmetry favors ferromagnetism. It is most use- 
ful to have a large DOS at the lower band edge in the hole 
doping picture to stabilize the Nagaoka state. Note that 
this is not equivalent to the well-known Stoner criterion 
Up(sf) > 1 which concerns only the DOS at the Fermi 



levels. The best situation is to have a strongly diverging 
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singularity at the lower band edge or close to it as we 
found in the investigation of the t-t' model with tunable 
DOS api- as-jwas also observed previously for /cc-type 
lattice&jM 

Out results concerning the t-t' model extend previous 
onesEJ since we treat all ratios of t and t' and all fillings 
ri. Hlubina et al. focussed on the Fermi levels at the van- 
Hove singularity. Thus the Stoner criterion is at the basis 
of their investigation albeit it goes technically beyond 
this mean-field criterion. 

Herrmann and Nolting used a two-pole method (SDA: 
self-consistent spectral density approach) to investigate 
ferromagnetism for the simple cubic, the dr^^oo hyper- 
cubic, the d = oo fee, and the bec latticalaO at zero 
and at finite temperature. Their qualitative findings for 
zero temperature are similar to ours. We like, however, 
to point out that the two-pole method they employ is in- 
deed a generalization of the Gutzwiller ansatz in (0) with 
/ = to finite temperatures and non-saturated magne- 
tizations. For T = and saturation it reduces to (^) 
with / = 0. Thus it is not astounding that the«-,found 
a good agreement to the results of Shastry et alr3. Our 
approaches go far beyond (^|) (barring the question of 
the extendability to finite temperatures). This can be 
seen for instance for the simple cubic lattice, where we 
found S CT = 0.237 well below <5 cr = 0.32BE1 Already 
Roth found by numerical calculation in Ihc variational 
subspace of RESO the number <5 cr = 0.24EjL.For the bee 
lattice one finds again that the SDA methooO reproduces 
the Gutzwiller result <5 cr = 0.324 for saturated ferromag- 
netism whereas RESO yields a considerably lower value 
of S CI — 0.239. Thus one is led to the conclusion that the 
SDA two-pole methocLcannot be exact as claimed in the 
strong coupling limitEl 
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APPENDIX A: THE FULL RESOLVENT ANSATZ 

RES3 

Computing the elements of the matrix luP — L (see 
(|l2|)) using Wick's theorem we obtain 

(Pi)k lk2 = n ■ S klk2 + \A\^ , (Ala) 
(Di)^ = [n(uj - e{k 2 )) , 

+ei(l + e(k 2 - k h )/(zt))} ■ S klk2 + |A|-V - ^b) (Alb) 
N ^i 2 = 

- e(k 2 ) - e[k x - k 2 - k h ) + s h )). (Ale) 

Since only the positions of the creation and the annihi- 
lation operator are interchanged between the states | 
and I'J'fc) one gets P2 and D2 from Pi and Di substitut- 
ing n by 6, lu by lu — U and e(k) by —s(k). As for Di 
( |l2| ) , D2 contains a diagonal matrix and a ^-independent 
part: 



D 2 = d 2 1 + (lu - U - e h ) uu 4 



(A2) 



with 



( d 2 1 )k 1 k 2 = S k,k 2 ■ [<K W -U + £{k 2 )) 

+ei(l - e(k 2 - k h )/(zt))] , (u) k = \A\-^ 2 (A3) 

for k, k x ,k 2 £ FS. Dj" 1 is known already from RESO 
( |l5| ) and we obtain the matrix elements of N + D^ 1 N for 
kl,k.2 G FS, q v q 2 S BZ \ FS as 

(N+D^N),^ = |A|- 2 £ (/( £l )5 £i£2 



lu - e h 



/(9 1 )/(£ 2 )' 



l + (w-Eb)ft(w) |A| J 
(e h - efe -g 1 - kh) + e(ki) - x 
(e b - e(k 2 -q 2 - k h ) + e(k 2 ) - s(g 2 )). 



(A4) 



A remarkable simplification occurs if terms like s(k — q) 
factorize to — e(k)e(q) / (zt) . This happens if, as for hy- 
percubic lattices, every component gives the same con- 
tribution to the sum over q. due to the symmetry of the 
Brillouin zone. In this case the corresponding matrix el- 
ement of the one particle Green's function is invariant 



under permutation of the components. Of course the 
bound state we are looking for has to display the same 
permutation symmetry. Making use of this argument and 
ass umi ng k b — 0, £b = —zt the product in the second line 
of (A4) simplifies to 



(ztf( 



zt 



1)( 



zt 



1)( 



f(£i) 

zt 



1)( 



zt 



1 



(A5) 



Carrying out the summation over q 1 and q 2 we obtain 

'e(ii) ,\ ( efe 



(N+D^N)^ 



|A|" 



zt 



zt 



1 j a 

(A6) 



with 



a 



7 _1 (e! + 2Szt - Sfl + (O + <) 2 G(f7) 

(ft + zi)[(ft + zf)G(ft)-<5] 2 - 
n + (Q + zt)G(Q) 



(A7) 



depending only on u) and ep, but not on the indices 
k x and fc 2 . We define the vector w by (w)k '■— 
|A|- 1 / 2 e(fc)/(z£) for k £ FS. Making use of this defi- 
nition, (|3l|), (A2), and (A6) the matrix B^ 1 reads 



b; 



d 2 1 + (lu + zt - U)uu + 
a(uu + — uw + — wu + + ww + ) 



(A8) 



The off-diagonal elements of B 2 1 do not depend explic- 
itly on k x and k 2 . In contrast to RESO (|12j), however, 
they are not overall constant but take specific values for 
each block of B 2 . To overcome this additional compli- 
cation we introduce the 2x2 matrix 

a 



A = 



ai 


03 




a 3 


0.2 





LU — Zt - 

—a 



U -a 
a 



(A9) 



and write B 2 1 as B 2 1 



2 1 - V + A y with j/+ := (u, w). 
In order to obtain B 2 we use an expansion trick similar 



to (14b) 



B 2 = d 2 {l-y+Ayd 2 )- 1 

= d 2 (l-y + Ayd 2 + y + Ayd 2 y + Ayd 2 -) 
= d 2 + d 2 y+A(l - BA)- X y d 2 , 

with B := y d 2 y + representing the 2x2 matrix 



B = 



bi 63 



u + d 2 u u + d 2 w 
u + d 2 w w + d 2 w 



(A10) 



(AH) 



While inverting the matrix 
1 -BA 



ci c 3 

C4 C 2 



1 - ai6i - a 3 b 3 
-ai6 3 - a 3 6 2 



-a 3 6i - a 2 6 3 
1 - a 2 6 2 - a 3 6 3 



(A12) 



19 



represents a simple algebraic task, the elements of B have 
to be computed by nume rical integration. In analogy to 
h(w) — jj+diu (see (16a)) we introduce 



h(oj) := u + d.2 u 

= lAr 1 J2 \^-U)+le{k) + e 1 ]- 1 (A13) 



fceFS 



with 7 := 5 — e±/(zt). Just as h{ui), h(uj) reduces to an 
integral over the DOS: 



h{uj) = 
with 

G(y) 



p(e) de 



S(lo — U) + ei + 7e 



p(e) de 
y + e 



n := 



5(u - U) + ei 

7 



(A15) 



The symmetry of the DOS with respect to e = allows to 
map G(y) to th e integral G(y) already defined in ( |16b[ ). 
Following (All), the elements of the matrix B are given 

by bi = 7" ^(H), b 2 =2_-\zt)- 2 {e 1 ~nTl + n 2 G(n)), 
and 6 3 = {jzty^n- nG(fi)). 

To find the energy of the bound state we ha ve to solve 
the equation (w + B2«)" 1 = 0. Starting from (A10) and 
writing u formally as u — y + e 1 we obtain 



u+B 2 u = e+y (d 2 + d 2 yA(l - BA)" 1 ^ d 2 ) y+e, 
= e+ (1 + BA(1 - BA)- 1 ) Be ± 
= e+(l-BA)- 1 Be 1 



(A16) 



and, using (A12), we finally obtain the equation 

C4C3 — cic 2 = (A17) 

for the lower edg e of t he spectrum of (u>). After in- 
serting all terms (A17) takes the form 



Pl ■ G(ft) - p 2 = 



(A18) 



with 



Pi = 0-7(0 + zt) 2 -(u + zt-U) (jizt) 2 - a(ei + nfl)) 
P2 = orf(n(Ti+2xt + ei) + n 2 (w + zt - {/)) + 7 2 (zt) 2 . 

Making use of the identities oj + zt — U = 7<5 _1 (r2 + zt) , 
ei = zt(n — 7) and introducing x '■= a l + 2 ^7i ( A18| ) 
simplifies to 



1 



5 + (SI + zt) G{Q) 



= 1 - 



2*x 



a(fi + zt) 



(A19) 



From the definition of a (A7) we derive the expression 



X = (fi + jrt) 1 



for We define in analogy 



(A20) 



X:= (ft + zt) 1 



1 



(5 + (SI + zt) G(Sl) 



(A21) 



and write ( A19| ) as x = ztx/ct. The eliminination of a 
finally leads to the simple result 



_L = 1 + 2 

zt x X 



(A22) 



f-'GCn), (A14) 



The Nagaoka instability line U CI (8) is obtained by as- 
suming lo = £p for a given Fermi energy, calculating 7, 
7, and x an d solve ( |A22|) n umerically with respect to 
X- Note that U enters ( |A22| ) solely via SI and hence via 
X- To compute the optimum spin flip energy for RES3 



for fixed values U and 5, we solve ( A22 ) with respect to 
to and subtract the Fermi energy ep from the solution 

m{u,8). 



APPENDIX B: GENERAL UNFRUSTRATED 
LATTICE 

In this appendix it will be shown t hat t he fo rmulae 
derived in section |l| and the formulae ( A20 - A22] ), apply 
to all unfrustrated, isotropic, homogeneous lattices with 
nearest neighbor hopping. In this context, 'homogeneous' 
means that all sites are equivalent; 'isotropic' means that 
all bonds in all directions are equivalent. 'Unfrustrated' 
means that the state cq := |A| -1 / 2 £\ is an eigen state 
of the kinetic Hamiltonian with eigen energy e\> = —z\t\ 
where t is the hopping element as in (|l|) and z is the 
coordination number. This requires t > 0, hence the 
absence of frustration. Note that the lattice does not 
need to be a Bravais lattice. The Bethe lattice, however, 
is not unfrustrated for z > 1 in the abave sense since 
its lower band edge is £b = —2y/z — l|t|c3 and not Sh — 
-^1*1- 

Let us denote by c+ the creation operators which di- 
agonalize the kinetic energy 



(Bl) 



and by a+ the site diagonal creation operators. The uni- 
tary transformation between these two bases has the ma- 
trix elements / aj 



(B2) 



which implies the expectation values with respect to the 
Nagaoka state \AP) 



<c QT aJ T ) = j for e a > e F 



(a+ T C/3 T ) = /^,- for e f3 < e F 



(B3a) 

1 p.j «ji ^ cjn (B3b) 
The homogeneity required implies that 

\f a j\ 2 S(uj - e a ) = constant (B4) 
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on the lattice, i.e. it does not depend on j. Unitarity 
yields furthermore 



(B5) 



First we address RESO with the ansatz (e a > ep) 

^a-^a^a+mf+j. (B6) 



The resulting matrix elements are obtaine d by W ick's 
theorem and re-expressed with the help of (B3a- 



Pa', a — ^^a',a ^ \fa',j\ \fct,j\ 3 



~fia',a t ^ ' |/a',i| |/q 

(Li) 



(B7a) 

(B7b) 
(B7c) 



The matrix inversion to be solved is 



(wP - L)" 1 = (D 



N 



D a > = <5 Q -, a (n(u - e„) + ei + (e 1 /zt)s a ) 1 



l (hi) 



1&) 

(B8b) 
Be) 



It can be re-expressed with the help of the matrices M, 
A, and V 



Mi 



E 



|/a,i| 2 |/a,j| 



n(uj - e a ) + ei + (ei/zt)e a 
V a ,j := |/a,j I 



(B9a) 
(B9b) 
(B9c) 



where the 5 are all spatial vectors connecting nearest 
neighbors. One obtains 



(uP - L)" 1 = D DVA E(-MA) n V+D 



\n=0 



(BIO) 



The key observation at this stage is that the vector u 
with uj = |A| -1 / 2 is an eigenvector both of the matrices 
M and A. The c orre sponding eigenvalue for M is found 
with the help of (Q 



|/a,i| |/q 



L A - J . n(uj — E a ) + ei + (ei/zt)e a 



(Bll) 



which simplifies due to ( ]B5| ) in the end to the form (16a). 
The correspondi ng e igenvalue of A is u> + zt = u> — Eb- 
So the series in (BIO) yields a vanishing denominator for 



= 1 + (uj — Eb)h(uj). Thus we derived ( fDj ) for a much 
broader class of lattices. 

The equations for RES1 and RES2 follow in analogy to 
the derivation in sect. |l|B. The ansatz RES1 is identical 
to (JlJ) for fc b = and the additional matrix elements are 
the same as in (^0|) once £k is replaced by e a . An impor- 
tant point to note is that the homogeneity (B4) ensures 
that N a couples indeed to the constant eigenvector u 



(V+DN)j = (S - h{uj)n{io - s h )) Uj 



(B12) 



for which the series summation in ( plO[) was achieved. 

For the ansatz RES2 we work with (|25|) for k h = and 
find the matrix elements ( p6| ) after replacing Ek by s a . 
Using 



(V+DN), = y Uj 



(B13) 



with y as in (]27cj), we obtain again ( J28| ) as condition for 
the variational spin flip energy. 

Let us now turn to RES3. We use {er < £f) 



*/3=E a iT c /3T a J t J- Ar ')^ : 



(B14) 



in analogy to (|29|) for the doubly occupied states. The 
matrices Di for $ Q and D2 for as in ( ^0|) are given 
by 

(Di) a >, a = 5 a >, a {n(uj - E a ) + a + e\e a /{zb)) + 

^Ei/«'.zi 2 i^i 2+i Ei/«Mi 2 i/^i 2 . ( B15a ) 

i (I'D 

(D 2 )/j',,3 = 8p>,p(5(u> - U + e ) + e x - eiep/(zt)) + 



'Ei/^ii 2 i/wi 2 +*Ei^.ii 2 i/« 



(B15b) 



where Di can be read off from (B8b,B8c) and D2 is 
analogous for the states $>p. 

The matrix N, which couples the doubly and the non- 
doubly occupied subspaces (see (30)), is obtained again 
via Wick's theorem and with (B3b,B4) 



dd) 

+*E(/a J </a,i/^/ W -|/a,i| a |//Ji| 2 ) • (B16) 

(id) 

In order to re-express the inverse matrix B2 = (D2 — 
N+Di^N)- 1 we define 



(Ci) 



(B17a) 
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P',0 

(B17b) 

%,£.';j,£ —tSi,j(3s',0 — <%o) + tSi _ 8> ,j 5g , _ 5' ( 1 — <5^o) 

- ty^J i+ s» ,j6s> ,063,0 , (B17c) 

5" 



Vj,5 ; /3 :— f/3,j+gfp,j 



(B17d) 



where the spatial vectors 5, 8 , and J" link nearest neigh- 
bors or equal 0. The result is 



B 2 = D 2 1 

OO 

D 2 " 1 V + E + CiE(C 2 E + CiE)"VD 2 " 1 (B18) 



n=0 



where we once again focus on the geometric series. In 
slight extension of the situation for RESO-2 we do not 
guess one common eigenvector of E, C±, and C 2 but 
a two-dimensional subspace spanned by u and v. The 
vectors are defined by m,o := |A| -1 / 2 and zero otherwise, 
and by v^gjto '■= (z|A|) -1 / 2 and zero otherwise. 
Straightforward calculation shows 



Eu = —ztu + \fztv_ 
~Ev = — \fztu + tv 

which can be summarized in 

E = —zt a b + 



(B19a) 
(B19b) 



(B20) 



with a:~u — v/ yjz and b := u + vj \fz. The matrix ele- 
ments of Ci with respect to u and v in obvious notation 
are 



' rv f rv 



1 + [w - s h )h 



a',a 

= — (h 2 - ^i^-gb) 
zt 2 \ 1 + (u) - £b)ho J 



c " = c ™ = 75fA|S I °' (Dl) "' 

1 Of ,Ot 

-1 hi 



(B21a) 



(B21b) 



x/zt l + (w — Eb)h 



(B21c) 



h n := 



e n p{e)de 
Vt- 1 e 



(B22) 



It is useful to keep the following relations in mind 



h\ = (—5 + (nui + ei)h Q )/-f 
hi = (ei + (nuj + ei)hi)/j . 



(B23a) 
(B23b) 



For C 2 very similar equations are derived after replacing 
a by (3, i.e. by changing the summation over the unoccu- 
pied levels to a summation over the occupied levels 



ho 



l + (u-U-s h )h ' 

c vv _ J_ (r _ h\(uj-U -£b) 

2 zt 2 \ 2 1 + (w - £7 - e b )h 

CUV g~\ vu 1 ^1 
2 — — ~~i=i~. ; 77 rr~ 

1 + (w — £/ - £b)fto 

where the generalization of ( A.13Q (ho = h) 



(B24a) 
(B24b) 
(B24c) 



e n p(e)de 



is used. The following relations hold 

hi = {n-{5{uj-U) + e 1 )h )h , 
h 2 = (ei - (S(uj -U) + ei)/n)/7 



(B25) 



(B26a) 
(B26b) 



Due to the particularly simple form of E in (B20) all 
we need to do is to calculate 



Cl 



C-2 



a + Cia = 



zt 



(ho 



2hi 
~zt 



l + e h )h 

+—(h 2 + (w - e h )(h 2 h - h 2 ))) , 
zt / 

zt f T 2hi 

~zT 



(B27a) 



b + C 2 b = 



l + (u-U- e b )h V 
+ ^- t (h + {lo-U- e h )(h 2 h - /i 2 ))) . (B27b) 

The condition for the singularity of ( Bl§| ) reads now 

1 = c lC2 (B28) 



which is equivalent to (A22) as can be shown by some 
tedious, but straightforward calculation. Thus we have 
completed the proof that the equations for RESO-3 de- 
rived in the main text for hypercubic lattices hold for all 
unfrustrated, isotropic, homogeneous lattices with near- 
est neighbor hopping. Only the coordination number and 
the DOS enter the evaluation of the RES ansatzes. 



where we use the generalization of (16a) (ho = h) 
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APPENDIX C: TRIANGULAR LATTICE 

For the triangular lattice with t < the lower band 
edge £b = — 3 j*| is reached at fc b = (47r/3, 0). Since 
kb 7^ Qj the integral 



h(cu) = 



i[uj - e(k)} + ei [1 - e{k - k h )/{zt)] 



fc£BZ\FS 

(CI) 

cannot be mapped onto a one-dimensional integral over 
the DOS but has to be evaluated explicitly in momentum 
space. 

The optimum spin flip energy for RESO for a given hole 
density 5 follows from the solution loq of the equation 
1 + (cj — Eb)h(w) = as Aeoo(<5) = ujq — £f (see sect. ||). 
To obtain the Fermi energy corresponding to the critical 
hole density S cr the equation 1 + (ep — £b)M £ F) = has 
to be solved numerically. 



For RES1, (23) holds also for fc b ^ 0, since 



\M 1 J2k £ (k — ^b) = £ bei/(zt) due to the symmetry 
of the lattice. Calculating N+U^N for RES2, however, 
the integrals 

e n (k) 



n[w - e(k)] + ei[l - s(k - k h )/(zt)} , ieBZ \ FS 

(C2) 



which for k h = simplify to ( B22 ) have to be computed 
for n = 1,2. Although the outline of the derivation re- 
mains unchanged, this causes some differences in the ana- 
lytic expressions for the optimum spin flip energy and the 
Nagaoka instability line compared with the case k b = 
(see sect. ||). 

Evaluating the full resolvent ansatz RES3, the product 
in the second line of (A4) can be written as 



4 




(C3) 



making use of e b = ~ z \t\/^- The permutation symmetry 
with respect to the primitive lattice vectors which is es- 
sential for the factorization e(k — q) = —e(k)e(q)/(zt) 
holds also for the triangular lattice. The matrix 
N+D^N is calculated to be 

N + D^ 1 N = aiuu + + a 2 (uw + +wu + ) + a 2 ww + (C4) 

with 



a-i = el - 2e h hi + h 2 = H[e h h - h x f 



(C5a) 



a 2 = e b h - -Sbhi 



1 



-H 



-2 1.2 



-Ebhhi 



-A 



"3 = £b h = + -TT- h 

I lb 



-H 



£ l h2 - IjZbhhx + —hj 



(C5b) 



(C5c) 



In ( |C5a| ) - ( C5c ), H is a short-hand notation for (w — 
£ b)/[l + (w — £b)h(u>)]. The method developed in ap- 
pendix ^ to calculate m + B 2M is applicable also for the 
triangular lattice up to eq. ( A17| ) which yields the opti- 
mum spin flip energy for RES3. 

The elements of the 2x2 matrices A and B are given 
by ai = «i - (oj - e h -U), a 2 = a s , a 3 = a 2 , &i = h{w), 
b 2 = /i 2 (w)/e b , b 3 = hi(ui)/sb with hi(uj) and h 2 (uj) 
defined in analogy to ( |C2| ) as integrals over the Fermi 
sphere. 



APPENDIX D: KAGOME LATTICE 

To prepare the derivation of the ansatzes RESO - RES3 
for the frustrated Kagome lattice we diagonalize the one- 
particle problem explicitly. 




FIG. 18. Segment of the kagome lattice. The vectors are 
used in the main text. The numbers refer to the three sites 
in each unit cell of this non-Bravais lattice. 

Since we deal with a non-Bravais lattice with three sites 
per unit cell we have to solve a 3x3 eigenvalue problem 
with 



faj = exp(ikj) 



(Dl) 



where r(j) g {1,2,3} denotes the sub-lattice to which 
site j belongs. The one-particle Hamiltonian acting on 
4> a T becomes 
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/ 



h(k) = 2t 





COS(=^ L ) 

\ cos(= i ) cos( 



cos(fcn 1 /2) cos(=§ 2 ) 

C0S (M^ 1 ) ) 







) 



o 



(D2) 



where we used the unit vectors 



and 
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as shown in 



fig. ([IS]). The secular equation of (D2) is 

= (-2f + A) (A 2 + 2tA - 2t 2 + te A (k)) , (D3) 

where e^(k) is the triangular dispersion J44|). From the 
secular equation one deduces ([l5|) easily. More important 
for the following is the observation that h(k) can be diag- 
onalized by an orthogonal, i.e. real, transformation since 
it is real symmetric. Thus the phase of f a j is completely 



given by the plane wave factor exp(ikj) in (Dl). 

Since we wish to treat the frustrated case (t < 0) 
we modify the ansatz (B6) by introducing an additional 
phase factor Xj depending only on the sub-lattice and 
being unity on sub-lattice 1, exp(27ri/3) on sub-lattice 2, 
and exp(— 27ri/3) on sub-lattice 3 



(D4) 



The resultin g matrix elements Pq,< jQ , and L-|- Q / iQ are the 
same as in ( B7a , B7b ) since the phase factor cancels at 
each site. But Lj Q / Q does change into 



J l a' ,a 



2zt ' 



,a — t ^ \ fa',i\ 2 \fa,j\ 2 Xf\ J 

(hi) 



(D5) 



The change in the second term is obvious. The change in 
the first term A\ is less trivial. In a first step one obtains 



Zl k,S ~ 



Transforming the terms of the sum like 5 
j — ► j + 5 leads to 



(D6) 
-S and 



(D7) 



= fajfaJ+lW+S ■ 

The last equality holds since <fi a ,T(j) m (Dl) is rea l- Hence 
only the real part of A^A - +s in (D6) matters. It is -1/2 



leading thus to the first term i n (|D5|). 

From the matrix elements ( |B7"a , B7fc|. D5 ) we find the 



relations which are analogous to ( B8b| ,B8c 



D9^|p9b| ) 



D Q ' : a = 6 a ',a(n(u - e a ) + e.\ - (e 1 /2zt)e a ) 1 , (D8a) 

i 

+ tJ2\U',i} 2 \f a , j \ 2 XtX j , (D8b) 

(id) 



\ \fa,i\ \fa,j\ 

M t i ■= > —, \ — ~ 7 — , n x , (D8c) 

h3 - ^ n{w - e a ) + e x - (ex/2zt)e a ' v ; 

A i-i : wS U ■ 'E'V^i'V -V, • (D8d) 



The vector u is again an eigenvector of the matri ces M 
and A. Its e igenvalue for M is in analogy to (Bll) iden- 
tical to (16a) with the adapted definition 



7K = n + ei/(2zt) 



(D9) 



The eigenvalue for A is u) — zt/2 = to — £b as before. 
So the series in ( B10D yields a vanishin g de nominato r fo r 
= 1 + (uj — eh)h(uj) with h(u>) as in (16a) with 7 (16c) 
replaced by 7k (D9). So the DOS, the lower band edge 
£b, and 7k are the only quantities to be changed in order 
that RESO (15 applies to the frustrated kagome lattice. 
For RES1 the ansatz reads 

:= |A|-V2^exp(ifcO)4%4l^'>A/- (D10) 



in extension of (jig). The resulting condition is identical 
to ( p3| ) with the adapted 7k in (DS) and, of course, £b = 
zt/2. 

For RES2 the ansatz reads 



|* a ) := lAr 1 ^ exp(zfc bi )a+a 3T a+|AA')A_+ (Dll) 

<ij> 



yielding again condition (|28|) with the adapted quantities, 
in particular 7 K := e\ + e<ili2z£). 

The ansatz for the doubly occupied states in RES3 is 
the extension of ( |B14[ ) 



(D12) 



The relations analogous to (B15a,B15b B16 |B17c ) read 
(Di)o',o = S a > ;a (n(uj - s a ) + ei - e 1 e a /(2zt)) + 

i (id) 
(D 2 )/3',/3 = V,/3 (^ - t/ + ep) + ei + e 1 e p /{2zt)) + 

^Efci 2 i/^ 2+t Ei/^ii 2 i/«i 2A i A i' ( D13b ) 



Na.,/3 = *E (/aj/a :i l//3:i| 2 ~ \f<x,i\ 2 ff3,j_f/3,i y 

id 

t E G^fiiLJtdfpA - \Ui?\fp£K\ ), (D13c) 



-t E ^i+5",j^£',0%,oA£ Aj 



(D13d) 
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So far the analogy to the treatment of unfrustrated 
lattices is perfect once the different form of £b, 7k, and 
of 7k = S + e i / (2zt) is taken into account. In particular 



the formulae (|B21 , B24 ) for the matrice s Ci a nd C2 carry 



over. But due to the differe nt fo rm of (D13d) the matrix 
E is changed compared to ( B19 ) 



Eu = zt/2u + y/ztv 
Ew = —s/ztu — t/2v 



(D14a) 
(D14b) 



Triangular lattice 

PA(e) = iy^tir 2 )- 1 ■ K [zx/zq] 
For t > 0, zq and Zx are given by 



(E4a) 



3 + 2^/3 - 7ft - {e/{2t)f for 2t < e < 3t 
4^3 - 7ft for - 6t < e < 2t ' 

(E4b) 



=> E = t 



Z/2 -y/i 

j~z -1/2 



(D14c) 



actin g on (u, v) . This matrix is no longer singular as was 
E in (|B20| ) . Thus we stay on th e 2 x 2 matrix level. The 
singularity condition based on (B18) is 



= dot 1 - C 2 E+CiE 



(D15) 



which can be evaluated easily. This concludes the deriva- 
tion for the RES3 ansatz on the frustrated kagome lattice. 



APPENDIX E: DOS FOR THE LATTICES 
CONSIDERED 

In this appendix we give the explicit formulae for the 
densities of states for the lattices discussed in sect. 111. 
K[m] stands £ar the complete elliptic integral of the first 
kind (see e.go). 



zx 



4^3 -e/t 



for 2t < e<3t 



3 + 2^3 - e/t - (e/(2t)) 2 for - 6t < e < 2t 



(E4c) 



For t < 0, the upper and lower intervals in ( |E4bD and 
(E4c) have to be replaced by — 3|t| < e < — 2\t\ and 
—2m < e < 6\t\, respectively. 



Honeycomb lattice 

Ph(e) = \e/t\ -p A {M-e 2 /t) (E5) 

Kagome lattice 

I S (e - 2t) + 1 11 + e/t\ ■ PA (3t - (e + t) 2 /t) 

(E6) 



Pk(e) 



Square lattice 



Pa (E) = (2\t\, 2 r l -K 



Psc{£) = 7T 1 



Simple cubic lattice 

112 du 



H VT^ 2 
ui = max(-l,-2-e/(2t)) 
u 2 = mm(l,2 - e/(2t)) 



p n (s + 2tu), 



Pbcc(e) 



bec lattice 

2 t m du 



r- ./| £ |/2 \/4u 2 - e 2 



pn(u), 



(El) 



(E2a) 

(E2b) 
(E2c) 



(E3) 



Hep lattice (t — —1) 

2 f 1 „ 

Phcp(e) = - dyE(y), 

k Jo 

with the integrand 

s(y) = 

p A (e--(2+ E -)j / 2 ) 



y/-2-e- 



■ e_ 



E + -(2+ £+ )y 2 -e_ 
p A (£-+(6-e-)^ 2 ) 

\/s+-e--(6-e-)y 2 



for e > 



and 



e± = e + 2 (1 ± y/F+Z) . 



(E7a) 



+V ^ + £+ ei -f 2 +^i^- e _ tor £ < U , 



(E7b) 



(E7c) 
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